It is often required to test a particular piece of equipment to determine whether it will properly function in an intended environment. One important factor is the vibrational energy generated in the intended environment. For example, the piece of equipment may be utilized in a space craft and therefore must survive the intense vibrational energy generated during lift off.
The actual vibrational energy experienced is the response to many independent vibrational sources and thus the distribution function of the amplitude over time is represented by a Gaussian distribution. Being Gaussian, this vibrational energy is completely specified by the amplitude distribution and the spectral power density function (PSD)(the first and second order statistics). In practice this reference PSD is often determined by recording the actual environmental vibrational energy and determining the PSD.
A well known system for subjecting a device under test (DUT) to vibrational energy is illustrated in FIG. 1. The DUT 10 is mounted on a shaker table 12 which receives a drive signal, y(t), from a driver 14. It would appear, at first glance, that the drive signal should be characterized by the reference PSD. However, the drive signal is applied to the shaker/DUT system which modifies the actual PSD of the energy applied to the DUT because of their combined transfer function. Accordingly, the drive PSD must be modified so that the PSD characterizing the actual motion of the sample is equal to the reference PSD.
A feedback system for generating a drive signal for subjecting the DUT to vibrational energy characterized by the reference PSD is depicted in the remainder of FIG. 1. A transducer 16 is coupled to the DUT 10 to convert the vibrations of the DUT 10 into an electrical signal. The output of the transducer is coupled to an analog to digital unit (ADC) 18 subsystem that may include signal conditioning such as filters, which has its output coupled to PSD estimation unit (S(.omega.)) 2020. An equalization unit 22 has one input coupled to the output of the PSD estimator (S(.omega.)) 20 and the second input connected to a reference spectrum, S(.omega.), and its output coupled to the input of a frequency domain randomization unit 24. An inverse Fourier transform unit (IFT) 26 has its input coupled to the output of the frequency domain randomization unit 24 and its output coupled to a window function unit 27. The output of the window function unit 27 is coupled to the input of a time domain randomization unit 28. The output of the time domain randomization unit 28 is coupled to the input of the driver 14 which includes a DAC, power amplifier, and may include additional signal conditioning. In practice, a given system may omit the time domain randomization unit 28, the frequency domain randomization unit 24, the window function unit 27, or some combination thereof.
The basic purpose of the Random Vibration Control System is to have the device under test (DUT) experience a prescribed vibration as sensed at one or more points on the DUT by an appropriate transducer (accelerometer, typically). Said prescribed vibration is usually required to be Gaussian (normal) which implies that the first and second order statistics are necessary and sufficient for complete specification. The first order statistic is simply the amplitude probability density function ##EQU1## where .sigma. is the RMS value of the random process. The second order statistic describes the frequency content or spectrum of the normally distributed signal which controls the correlation (or covariance) between successive values of the random signal. Either, the PSD, S(.omega.), or its Fourier transform, the auto correlation function, .phi.(t), are usually used to specify this statistic because these two functions constitute a Fourier transform pair EQU S(.omega.) .rarw..fwdarw..phi.(t)
Because digital signal processing techniques are used in these control systems, uniformly sampled time domain and frequency domain functions are used in the control loop: analog to digital converters and discrete Fourier transform usually providing the means for obtaining sampled functions.
It is important to recognize that an essential feature of these normal random signals is that the discrete Fourier transform X(k) of a sample segment of x(n), consisting of N samples duration given by ##EQU2## is a complex function of the frequency index, k, such that EQU X(k)=P.sub.x (k)+jQ.sub.x (k)
where N-1 ##EQU3## Note that if the set of values {x(n)} are randomly distributed, the real part P.sub.x (k) and the imaginary part, Q.sub.x (k), being weighted sums of these randomly distributed samples must also be random. Further, if x(n) is normally distributed, P.sub.x (k) and C.sub.x (k) will be, for certain, normally distributed. Thus, the Fourier transform of a normally distributed signal is a complex random function with normally distributed real and imaginary parts. P.sub.x (k) and Q.sub.x (k) are, also, uncorrelated because the weights in cos 2.pi.kn/N and sin 2.pi.kn/N used to obtain each part are orthogonal (uncorrelated).
If the complex function, X(k), is expressed in polar coordinate form, EQU X(k)=.vertline.X(k).vertline.e.sup.j.theta.(k)
where ##EQU4## then .theta.(k) can be shown to be uniformly distributed on the interval 0.ltoreq..theta..ltoreq.2.pi., while .vertline.X(k).vertline. is Rayleigh distributed.
In order to establish control of the vibrations, the spectrum of the vibrations must be estimated and compared with the desired spectrum so that the spectral content of the drive signal may be appropriately tailored to produce the desired result. Conceptually, it is a straightforward problem: estimate the transfer function by estimating the transfer function between the input to the DAC 14 and the output of the ADC 18 from which the necessary drive signal may be determined. In practice, problems may arise because the vibrations are random, the system may be noisy and may also be nonlinear.
Because the sensed response, z(n) is random, the Fourier transform, in accordance with the previous discussion, is also random. This means that each of the individual estimated spectrums or periodograms will be a random function. Specifically, the spectral estimator of the response signal that we desire is ##EQU5## where ##EQU6## or EQU Z.sub.m (k)=P.sub.zm (k)+jQ.sub.zm (k),
.pi. is the sampling interval and z.sub.m (n) is the m.sup.th response data frame. Hence ##EQU7## where p.sub.zm.sup.2 (k) and Q.sub.zm.sup.2 (k) are each squared Gaussian variables that are, therefore, chi-squared distributed with one degree of freedom each. The sum of the terms ##EQU8## known as the "raw" spectral estimates or periodograms, is also chi-squared distributed but with two degrees of freedom, and has a relative variance of unity with respect to its true mean value. The purpose for averaging over an infinite ensemble of such terms is to estimate S.sub.z (k) by reducing these variances. In practice a finite weighted sum of "raw" spectral estimates, S.sub.zm (k), are used as a basis of control. For example EQU S.sub.z,m+l (k)+.alpha.S.sub.z,m (k)+(l-.alpha.) S.sub.z,m+l (k)
where
S.sub.z,m (k) is the m.sup.th spectral estimate; ##EQU9## .alpha. is a "discount" averaging factor (0 &lt;.alpha.&lt;1). Thus, this spectral estimate will have errors (variance) because of the finite number of degrees of freedom used in the estimate. The variance decreases inversely with m. Additional variance may be introduced by such non-linear PSD estimating techniques as windowing of the sampled response data, z(n), for purposes of reducing leakage in the spectral estimate. PA0 1. Dynamic range is limited because of the windowing constraint. PA0 2. True phase and frequency randomization are lacking. PA0 3. Even under the best circumstances, the second order statistics are non-stationary because the overlap and windowing operations result in a correlation function that is a function of two times, .phi.(t.sub.1,t.sub.2), rather than the stationary function which is only dependent on the time difference T=t.sub.1 -t.sub.2, or .phi.(t.sub.1 -t.sub.2)=.phi.(T).
Because the system may also be noisy and non-linear, a stochastic approximation algorithm is used that introduces less than the total correction implied by the latest average response spectrum. Thus, a typical correction for the drive Fourier transform magnitude may take on the form ##EQU10## where 0&lt;.epsilon.&lt;1/2 is the stochastic weighting factor. This procedure results in a deterministic magnitude function whereas what we desire is EQU X.sub.m+1 (k)=.vertline.X.sub.m+1 (k).vertline.e.sup.j.sup..theta. m+l.sup.(k) =P.sub.x,m+1 (k)+jQ.sub.x,m+l (k)
where .vertline.X.sub.m+l (k).vertline. is Rayleigh distributed, .theta..sub.m+1 (k) is uniformly distributed or equivalently P.sub.x,m+1 and Q.sub.x,m+1 are independently normally distributed in order to produce a proper Gaussian distributed drive function, x.sub.m+1 (n) after taking the inverse Fourier transform of X.sub.m+1 (k).
Prior art has established that an acceptable approximation to the normal drive signal may be achieved by only randomizing X(k) by introducing a random phase function .theta.(k). See U.S. Pat. No. 3,710,082. Because the central limit theorem asserts that the sum of a large number of random variables tends to be Gaussian, almost independent of their particular distributions, the inverse Fourier transform, x.sub.m+1 (n), will also be approximately Gaussian. This procedure has been dubbed "frequency domain" randomization even though only the phase is randomized and not the amplitude also.
A further variant was also introduced by Sloane et al., see U.S. Pat. No. 3,848,115 that has been popularly referred to as time domain randomization that is based on randomizing phase only by random circular rotation of x(n). This procedure also only randomizes phase, and not both phase and amplitude.
The time domain function, x.sub.m (n), is a finite length array of N data points that could be concatenated with other such arrays to produce a continuous stream, x(m,n) such that ##EQU11## where 0&lt;n&lt;N-1. However, it was found in Sloane et al. that the uniform weighting implied by construction of x(m,n) introduced unacceptable spectral leakage that limited the ability to control systems exhibiting large sharp resonances or being able to control when the desired reference spectrum, S(k), had large dynamic range characteristics. As a result, non-uniform windowing of x.sub.m (n) was introduced in order to control this spectrum.
FIG. 2A and 2B depict systems utilizing phase domain and time domain randomization, respectively, followed by windowing and overlap and add.
From a practical point of view, non-uniform windowing by itself could not be arbitrarily applied to each segment, x.sub.m (n), because concatenation of such a sequence, ##EQU12## where w(n) is the window function, would produce a non-stationary output sequence whose instantaneous power was a periodic function proportional to s.sup.2 (n) in violation of the usual requirement that, at steady-state or at convergence of control to the proper spectrum, the drive function must also be a stationary normal process. This led to the introduction of overlapping in addition to windowing. A further implicit constraint was thus established on the choice of window functions because, for example, a 50% overlap of windowed arrays required that where c is a constant. For the 50% overlap case, this restricted the choice of window to the half-sine window ##EQU13## Unfortunately, this constraint on the window function restricted the choice to one of the least effective of all window function. Thus, we conclude that the state of the art is deficient on several